rm(list=ls());gc()
##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 428038 22.9     898738   48   648454 34.7
## Vcells 771478  5.9    8388608   64  1661428 12.7
load("G:/Mi unidad/paperestallido/Definitive_models_2021.RData")
#xaringan::inf_mr()
#setwd("C:/Users/CISS Fondecyt/Pictures")
#C:\Users\CISS Fondecyt\Pictures
#arriba puse algunas opciones para que por defecto escondiera el código
#también cargue algunos estilo .css para que el texto me apareciera justificado, entre otras cosas.
local({r <- getOption("repos")
       r["CRAN"] <- "http://cran.r-project.org" 
       options(repos=r)
})

`%>%` <- magrittr::`%>%`
copy_names <- function(x,row.names=FALSE,col.names=TRUE,dec=",",...) {
  if(class(ungroup(x))[1]=="tbl_df"){
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)    
        }
  } else {
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)       
  }
 }
}  

unlink('*_cache', recursive = TRUE)
if(!require(pacman)){install.packages("pacman")}

pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes
knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
#dejo los paquetes estadísticos que voy a utilizar

if(!require(plotly)){install.packages("plotly")}
if(!require(lubridate)){install.packages("lubridate")}
if(!require(htmlwidgets)){install.packages("htmlwidgets")}
if(!require(tidyverse)){install.packages("tidyverse")}
if(!require(gganimate)){install.packages("gganimate")}
if(!require(readr)){install.packages("readr")}
if(!require(stringr)){install.packages("stringr")}
if(!require(data.table)){install.packages("data.table")}
if(!require(DT)){install.packages("DT")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(lattice)){install.packages("lattice")}
if(!require(forecast)){install.packages("forecast")}
if(!require(zoo)){install.packages("zoo")}
if(!require(panelView)){install.packages("panelView")}
if(!require(janitor)){install.packages("janitor")}
if(!require(rjson)){install.packages("rjson")}
if(!require(estimatr)){install.packages("estimatr")} 
if(!require(CausalImpact)){install.packages("CausalImpact")}
if(!require(textreg)){install.packages("textreg")}
if(!require(sjPlot)){install.packages("sjPlot")}
if(!require(foreign)){install.packages("foreign")}
if(!require(tsModel)){install.packages("tsModel")}
if(!require(lmtest)){install.packages("lmtest")}
if(!require(Epi)){install.packages("Epi")}
if(!require(splines)){install.packages("splines")}
if(!require(vcd)){install.packages("vcd")}
if(!require(astsa)){install.packages("astsa")}
if(!require(forecast)){install.packages("forecast")}
if(!require(MASS)){install.packages("MASS")}
if(!require(ggsci)){install.packages("ggsci")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(dplyr)){install.packages("dplyr")}
if(!require(ggforce)){install.packages("ggforce")}
if(!require(imputeTS)){install.packages("imputeTS")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(SCtools)){install.packages("SCtools")}
if(!require(MSCMT)){install.packages("MSCMT")}
# Calculate the number of cores
no_cores <- detectCores() - 1
cl<-makeCluster(no_cores)
registerDoParallel(cl)

Sys.setlocale(category = "LC_ALL", locale = "english")
## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

0.1 Final Plots

Black lines are the observed trend for each outcome, red lines are the estimated trends through Bayesis structural times-series model, blue areas are the 95% credible interval from estimates, and vertical line is the onset of social protests on October 18th, 2019


library(cowplot)

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-1.2 #.8
angle_x<-45
Sys.setlocale(category = "LC_ALL", locale = "english")
## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
hosp_trauma_plot_red<-
plot(impact3d_hosp_trauma, "original")$data %>% 
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
    ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,125))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
        theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = ggtext::element_markdown(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = ggtext::element_markdown(size=size_title),
          plot.caption = ggtext::element_markdown(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot_red<-
plot(impact3d1_hosp_resp, "original")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,125))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = ggtext::element_markdown(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = ggtext::element_markdown(size=size_title),
          plot.caption = ggtext::element_markdown(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

cons_trauma_plot_red<-
plot(impact3d1_cons_trauma, "original")$data %>% 
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,1600))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot_red<-
plot(impact3d_cons_resp, "original")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,1600))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

rate_trauma_plot_red<-
plot(impact3d_ratio, "original")$data %>% 
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,400))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot_red<-
plot(impact3d_ratio_resp, "original")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,400))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
plot1b_red<-
plot_grid(
  cons_trauma_plot_red+theme(axis.text.x = element_blank()),
  cons_resp_plot_red+theme(axis.text.x = element_blank()),
  hosp_trauma_plot_red+theme(axis.text.x = element_blank()),
  hosp_resp_plot_red+theme(axis.text.x = element_blank()),
  rate_trauma_plot_red,
  rate_resp_plot_red,#,
  nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25)
)

ggdraw(plot1b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))
Appendix Figure 1. Trends of emergency department consultations and hospitalizations (2015-2019)

Appendix Figure 1. Trends of emergency department consultations and hospitalizations (2015-2019)

  jpeg("_figs/_FigS1_final.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot1b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2
  #add_sub(, "Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nBlue area= Prediction intervals.", vpadding=grid::unit(0, "lines"), y = .55, x = 0.03, hjust = 0,size=9)
#plot2 <- cowplot::ggdraw(grid.arrange(p14, p21,p212,p32,p42,p52,p57, ncol = 2, nrow = 4)) + 
  # same plot.background should be in the theme of p1 and p2 as mentioned above
#  theme(plot.background = element_rect(fill=NA, color = NA))

cairo_ps("_figs/_FigS1_final.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      ggdraw(plot1b_red)+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_FigS1_final.ps", plot =
          ggdraw(plot1b_red)+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)


Blue dashed lines are the estimated difference between predicted (the counterfactual) and observed trend in the ten-weeks before and after the onset of the social protests. Blue areas are the 95% credible intervals. Vertical dashed line signals the onset of social protests on October 18th, 2019.


#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8
angle_x<-45

hosp_trauma_plot2_red<-
plot(impact3d_hosp_trauma, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
  theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

hosp_resp_plot2_red<-
plot(impact3d1_hosp_resp, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_trauma_plot2_red<-
plot(impact3d1_cons_trauma, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_resp_plot2_red<-
plot(impact3d_cons_resp, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_trauma_plot2_red<-
plot(impact3d_ratio, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_resp_plot2_red<-
plot(impact3d_ratio_resp, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #añadí el 2021-05-11
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
  #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

vector_breaks_plot2b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")
## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
vector_dates_plot2b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))

plot2b_red<-
plot_grid(
  cons_trauma_plot2_red+scale_y_continuous(breaks=seq(-800,250,200), limits = c(-800,250)),
  cons_resp_plot2_red+scale_y_continuous(breaks=seq(-800,250,200), limits = c(-800,250)),
  hosp_trauma_plot2_red+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2_red+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2_red+scale_y_continuous(breaks=seq(-200,400,200), limits = c(-200,400)),
  rate_resp_plot2_red+scale_y_continuous(breaks=seq(-200,400,200), limits = c(-200,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(.95, .95), rel_heights = c(rep(1,2),1.3), scale = c(.95, .95, .95,.95,.95,.95)
)+
  draw_label("Weekly difference change in consultations/hospitalizations",x=0.005, y=.5, angle=90, size = 17)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

    ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
Figure 1. Weekly differences between predicted and observed outcomes in the 10 weeks pre and post-exposure periods

Figure 1. Weekly differences between predicted and observed outcomes in the 10 weeks pre and post-exposure periods

  jpeg("_figs/_Fig1_final.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2
cairo_ps("_figs/_Fig1_final.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      ggdraw(plot2b_red)+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_Fig1_final.ps", plot =
          ggdraw(plot2b_red)+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)


Blue dashed lines are the estimated cumulative difference between predicted (the counterfactual) and observed trend in the ten-weeks before and after the onset of the social protests. Blue areas are the 95% credible intervals. Vertical dashed line signals the onset of social protests on October 18th, 2019.


#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8
angle_x <- 45

hosp_trauma_plot3_red<-
plot(impact3d_hosp_trauma, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
    dplyr::filter(date>="2019-08-12") %>% 
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

hosp_resp_plot3_red<-
plot(impact3d1_hosp_resp, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_trauma_plot3_red<-
plot(impact3d1_cons_trauma, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_resp_plot3_red<-
plot(impact3d_cons_resp, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_trauma_plot3_red<-
plot(impact3d_ratio, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_resp_plot3_red<-
plot(impact3d_ratio_resp, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

Sys.setlocale(category = "LC_ALL", locale = "english")
## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
plot3b_red<-
plot_grid(
  cons_trauma_plot3_red+scale_y_continuous(breaks=seq(-4000,1500,1000), limits = c(-4000,1500)),
  cons_resp_plot3_red+scale_y_continuous(breaks=seq(-4000,1500,1000), limits = c(-4000,1500)),
  hosp_trauma_plot3_red+scale_y_continuous(breaks=seq(-200,200,100), limits = c(-200,200)),
  hosp_resp_plot3_red+scale_y_continuous(breaks=seq(-200,200,100), limits = c(-200,200)),
  rate_trauma_plot3_red+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3_red+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25), scale = c(.95, .95, .95,.95,.95,.95))+
 draw_label("Weekly difference change in consultations/hospitalizations",x=0.005, y=.5, angle=90, size = 17)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

ggdraw(plot3b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))
Figure 3. Cumulative difference between predicted and observed outcomes in the 10 weeks pre and post-exposure period

Figure 3. Cumulative difference between predicted and observed outcomes in the 10 weeks pre and post-exposure period

  jpeg("_figs/_FigS2_final.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot3b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2
cairo_ps("_figs/_Fig2_final.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      ggdraw(plot3b_red)+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_FigS2_final.ps", plot =
          ggdraw(plot3b_red)+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)

0.2 Consolidated Plot

# dejar una sola figura que contenga 6 líneas (una para cada outcome) con el cambio porcentual (para que todos los outcomes estén en la misma unidad de medida).

#_#_#_#_#_# CONS TRAUMA #_#_#_#_#_#

releff_cons_trauma<-
cbind(plot(impact3d1_cons_trauma, "original")$data,
      data15a64_rn[,c("date","did")])%>%
      dplyr::mutate_at(.vars=vars(response,mean,lower,upper), 
                   .funs = funs(`zw`=scale(.))) %>% 
      dplyr::mutate(diff=response-mean) %>% 
      dplyr::mutate(rel_diff=diff/mean) %>%   
      dplyr::mutate(diff_zw=response_zw-mean_zw) %>% 
      #ggplot(aes(y=rel_diff, x=time))+
     #   geom_line()+
     # geom_vline(xintercept = 253, color="red")+
     # theme_bw()
      dplyr::select(date,time,rel_diff,diff_zw) #%>% 
      #tail(20)
  #ggplot(aes(y=rel_diff, x=date))+
  #geom_line()
# impact3d1_cons_trauma
#  Relative effect (s.d.)   -14% (13%)    -14% (13%)   
#  95% CI                   [-40%, 12%]   [-40%, 12%]    

#_#_#_#_#_# CONS RESP #_#_#_#_#_#

releff_cons_resp<-
cbind(plot(impact3d_cons_resp, "original")$data,
      data15a64_rn[,c("date","did")])%>% 
      dplyr::mutate_at(.vars=vars(response,mean,lower,upper), 
                   .funs = funs(`zw`=scale(.))) %>% 
      dplyr::mutate(diff=response-mean) %>% 
      dplyr::mutate(rel_diff=diff/mean) %>%   
      dplyr::mutate(diff_zw=response_zw-mean_zw) %>% 
      dplyr::select(date,time,rel_diff,diff_zw) #%>% 
      #tail(20)
  #ggplot(aes(y=rel_diff, x=date))+
  #geom_line()
# impact3d_cons_resp
# Relative effect (s.d.)   -30% (30%)    -30% (30%)  
# 95% CI                   [-89%, 30%]   [-89%, 30%]    


#_#_#_#_#_# HOSP TRAUMA #_#_#_#_#_#

releff_hosp_trauma<-
cbind(plot(impact3d_hosp_trauma, "original")$data,
      data15a64_rn[,c("date","did")])%>% 
      dplyr::mutate_at(.vars=vars(response,mean,lower,upper), 
                   .funs = funs(`zw`=scale(.))) %>% 
      dplyr::mutate(diff=response-mean) %>% 
      dplyr::mutate(rel_diff=diff/mean) %>%   
      dplyr::mutate(diff_zw=response_zw-mean_zw) %>% 
      dplyr::select(date,time,rel_diff,diff_zw) #%>% 
      #tail(20)
  #ggplot(aes(y=rel_diff, x=date))+
  #geom_line()
# impact3d_cons_resp
# Relative effect (s.d.)   15% (5.6%)   15% (5.6%) 
# 95% CI                   [4%, 26%]    [4%, 26%]


#_#_#_#_#_# HOSP RESP #_#_#_#_#_#

releff_hosp_resp<-
cbind(plot(impact3d1_hosp_resp, "original")$data,
      data15a64_rn[,c("date","did")])%>% 
      dplyr::mutate_at(.vars=vars(response,mean,lower,upper), 
                   .funs = funs(`zw`=scale(.))) %>% 
      dplyr::mutate(diff=response-mean) %>% 
      dplyr::mutate(rel_diff=diff/mean) %>%   
      dplyr::mutate(diff_zw=response_zw-mean_zw) %>% 
      dplyr::select(date,time,rel_diff,diff_zw) #%>% 
      #tail(20)
  #ggplot(aes(y=rel_diff, x=date))+
  #geom_line()
# impact3d_cons_resp
#  Relative effect (s.d.)   -6.8% (19%)   -6.8% (19%)  
# 95% CI                   [-44%, 31%]   [-44%, 31%]     


#_#_#_#_#_# RATIO TRAUMA #_#_#_#_#_#

releff_ratio_trauma<-
cbind(plot(impact3d_ratio, "original")$data,
      data15a64_rn[,c("date","did")])%>% 
      dplyr::mutate_at(.vars=vars(response,mean,lower,upper), 
                   .funs = funs(`zw`=scale(.))) %>% 
      dplyr::mutate(diff=response-mean) %>% 
      dplyr::mutate(rel_diff=diff/mean) %>%   
      dplyr::mutate(diff_zw=response_zw-mean_zw) %>% 
      dplyr::select(date,time,rel_diff,diff_zw) #%>% 
      #tail(20)
  #ggplot(aes(y=rel_diff, x=date))+
  #geom_line()
# impact3d_cons_resp
# Relative effect (s.d.)   40% (14%)    40% (14%)  
# 95% CI                   [13%, 68%]   [13%, 68%] 


#_#_#_#_#_# RATIO RESP #_#_#_#_#_#

releff_ratio_resp<-
cbind(plot(impact3d_ratio_resp, "original")$data,
      data15a64_rn[,c("date","did")])%>% 
      dplyr::mutate_at(.vars=vars(response,mean,lower,upper), 
                   .funs = funs(`zw`=scale(.))) %>% 
      dplyr::mutate(diff=response-mean) %>% 
      dplyr::mutate(rel_diff=diff/mean) %>%   
      dplyr::mutate(diff_zw=response_zw-mean_zw) %>% 
      dplyr::select(date,time,rel_diff,diff_zw) #%>% 
      #tail(20)
  #ggplot(aes(y=rel_diff, x=date))+
  #geom_line()
# impact3d_cons_resp
#  Relative effect (s.d.)   59% (15%)    59% (15%)   
#  95% CI                   [29%, 88%]   [29%, 88%]
#horizontal
height_y<-18
height_x<-18
size_title<-18
line_size<-1.3
angle_x <- 45

Sys.setlocale(category = "LC_ALL", locale = "english")
## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
manualcolors<-c('indianred1','cornflowerblue', 'gray50', 'darkolivegreen4', 'slateblue2', 'firebrick4', 'goldenrod4')

releff<-
releff_cons_trauma[,c("date","rel_diff")] %>% 
  dplyr::left_join(releff_cons_resp[,c("date","rel_diff")],by="date") %>%
  dplyr::left_join(releff_hosp_trauma[,c("date","rel_diff")],by="date") %>%
  dplyr::left_join(releff_hosp_resp[,c("date","rel_diff")],by="date") %>%
  dplyr::left_join(releff_ratio_trauma[,c("date","rel_diff")],by="date") %>%
  dplyr::left_join(releff_ratio_resp[,c("date","rel_diff")],by="date")

colnames(releff)<-c("date","Cons_Trauma","Cons_Resp","Hosp_Trauma","Hosp_Resp","Ratio_Trauma","Ratio_Resp")

releff_melt<-
releff %>% 
  melt(id.vars=c("date")) %>% 
  dplyr::mutate(variable=
          dplyr::case_when(variable=="Cons_Trauma"~"Trauma Consultations",
                           variable=="Cons_Resp"~"Respiratory Consultations",
                           variable=="Hosp_Trauma"~"Trauma Hospitalizations",
                           variable=="Hosp_Resp"~"Respiratory Hospitalizations",
                           variable=="Ratio_Trauma"~"Trauma Hospitalizations per 1,000 consultations",
                           variable=="Ratio_Resp"~"Respiratory Hospitalizations per 1,000 consultations")) %>% 
  dplyr::filter(date>="2019-08-05")

releff_plot<-
releff_melt%>% 
  dplyr::mutate(variable=factor(variable)) %>% 
  ggplot(aes(x=date, y=value, color=variable))+ #, shape=variable
  geom_line(size=line_size)+
  #scale_shape_manual(name= "Outcome", values=c(1:6))+
  #scale_linetype_manual(name= "",
  #                      values= c("twodash","solid","dotted","dashed","twodash","solid"))+
  scale_color_manual(name= "", values=manualcolors)+
  theme_sjplot()+
  theme(legend.position = c(.2,.8),
        legend.text= element_text(size=18),
        legend.key.width = unit(4,"cm"),
        legend.key.height = unit(1.4,"cm"),
        legend.background = element_rect(fill = "white", color = "black"),
        legend.title=element_blank(),
        axis.text.y= element_text(size=height_y),
        axis.title.y= element_text(size=height_y),
        axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x))+
  labs(y="Percent change", x="")+
  scale_x_date(breaks = scales::date_breaks("4 weeks"), labels = scales::date_format("%d-%b"))+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  geom_vline(xintercept = as.Date("2019-10-21"), col = "darkred", lty = 2, size=1)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  guides(color = guide_legend(override.aes = list(size = 3) ) )

releff_plot
Figure 4. Relative Effects

Figure 4. Relative Effects

cairo_ps("_figs/_Fig1_final_AJPH.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      releff_plot+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_Fig1_final_AJPH.ps", plot =
          releff_plot+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)
 
jpeg("_figs/_Fig1_AJPH.jpg", height = 14, width = 23, units = 'in', res = 600)
  releff_plot+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2

0.3 ESP

library(cowplot)

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-1.2 #.8
angle_x<-45
Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
hosp_trauma_plot_red_sp<-
plot(impact3d_hosp_trauma, "original")$data %>% 
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
    ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,125))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
        theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = ggtext::element_markdown(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = ggtext::element_markdown(size=size_title),
          plot.caption = ggtext::element_markdown(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Hospitalizaciones por Causa: Traumatismos")

hosp_resp_plot_red_sp<-
plot(impact3d1_hosp_resp, "original")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,125))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = ggtext::element_markdown(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = ggtext::element_markdown(size=size_title),
          plot.caption = ggtext::element_markdown(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Hospitalizaciones por Causa: Sistema Respiratorio")

cons_trauma_plot_red_sp<-
plot(impact3d1_cons_trauma, "original")$data %>% 
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,1600))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Consultas por Causa: Traumatismos")

cons_resp_plot_red_sp<-
plot(impact3d_cons_resp, "original")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,1600))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Consultas por Causa: Sistema Respiratorio")

rate_trauma_plot_red_sp<-
plot(impact3d_ratio, "original")$data %>% 
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,400))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Hospitalizaciones por Causa: Traumatismos, cada 1.000 consultas")

rate_resp_plot_red_sp<-
plot(impact3d_ratio_resp, "original")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_date(date_breaks = "3 months",date_labels="%b %y",limits=c(as.Date("2015-01-05"),as.Date("2019-12-23")))+
  ylim(c(0,400))+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Hospitalizaciones por Causa: Sistema Respiratorio, cada 1.000 consultas")

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:


Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
plot1b_red_sp<-
plot_grid(
  cons_trauma_plot_red_sp+theme(axis.text.x = element_blank()),
  cons_resp_plot_red_sp+theme(axis.text.x = element_blank()),
  hosp_trauma_plot_red_sp+theme(axis.text.x = element_blank()),
  hosp_resp_plot_red_sp+theme(axis.text.x = element_blank()),
  rate_trauma_plot_red_sp,
  rate_resp_plot_red_sp,#,
  nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25)
)

ggdraw(plot1b_red_sp)+
   theme(plot.background = element_rect(fill=NA, color = NA))
Figura Suplementaria 1. Tendencias de consultas y hospitalizaciones por emergencias (2015-2019)

Figura Suplementaria 1. Tendencias de consultas y hospitalizaciones por emergencias (2015-2019)

  jpeg("_figs/_FigS1_final_ESP.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot1b_red_sp)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2
  #add_sub(, "Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nBlue area= Prediction intervals.", vpadding=grid::unit(0, "lines"), y = .55, x = 0.03, hjust = 0,size=9)
#plot2 <- cowplot::ggdraw(grid.arrange(p14, p21,p212,p32,p42,p52,p57, ncol = 2, nrow = 4)) + 
  # same plot.background should be in the theme of p1 and p2 as mentioned above
#  theme(plot.background = element_rect(fill=NA, color = NA))

cairo_ps("_figs/_FigS1_final_ESP.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      ggdraw(plot1b_red_sp)+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_FigS1_final_ESP.ps", plot =
          ggdraw(plot1b_red_sp)+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)
#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8
angle_x<-45

Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
hosp_trauma_plot2_red<-
plot(impact3d_hosp_trauma, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
  theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Hospitalizaciones por Causa: Traumatismos")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

hosp_resp_plot2_red<-
plot(impact3d1_hosp_resp, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Hospitalizaciones por Causa: Sistema Respiratorio")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_trauma_plot2_red<-
plot(impact3d1_cons_trauma, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Consultas por Causa: Traumatismos")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_resp_plot2_red<-
plot(impact3d_cons_resp, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Consultas por Causa: Sistema Respiratorio")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_trauma_plot2_red<-
plot(impact3d_ratio, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Hospitalizaciones por Causa: Traumatismos, cada 1.000 consultas")+
    #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_resp_plot2_red<-
plot(impact3d_ratio_resp, "pointwise")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  #añadí el 2021-05-11
  #2021-05-11
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Hospitalizaciones por Causa: Sistema Respiratorio, cada 1.000 consultas")+
  #Change in 2021-05-11
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

vector_breaks_plot2b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
vector_dates_plot2b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))

Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
plot2b_red<-
plot_grid(
  cons_trauma_plot2_red+scale_y_continuous(breaks=seq(-800,250,200), limits = c(-800,250)),
  cons_resp_plot2_red+scale_y_continuous(breaks=seq(-800,250,200), limits = c(-800,250)),
  hosp_trauma_plot2_red+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2_red+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2_red+scale_y_continuous(breaks=seq(-200,400,200), limits = c(-200,400)),
  rate_resp_plot2_red+scale_y_continuous(breaks=seq(-200,400,200), limits = c(-200,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(.95, .95), rel_heights = c(rep(1,2),1.3), scale = c(.95, .95, .95,.95,.95,.95)
)+
  draw_label("Diferencias semanales en consultas/hospitalizaciones",x=0.005, y=.5, angle=90, size = 17)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

    ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
Figura 1. Diferencias semanales entre resultados predichos y observados en las 10 semanas anteriores y posteriores a la exposición

Figura 1. Diferencias semanales entre resultados predichos y observados en las 10 semanas anteriores y posteriores a la exposición

  jpeg("_figs/_Fig1_final_ESP.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2
cairo_ps("_figs/_Fig1_final_ESP.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      ggdraw(plot2b_red)+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_Fig1_final_ESP.ps", plot =
          ggdraw(plot2b_red)+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)
#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8
angle_x <- 45

Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
hosp_trauma_plot3_red<-
plot(impact3d_hosp_trauma, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
    dplyr::filter(date>="2019-08-12") %>% 
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Hospitalizaciones por Causa: Traumatismos")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

hosp_resp_plot3_red<-
plot(impact3d1_hosp_resp, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Hospitalizaciones por Causa: Sistema Respiratorio")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_trauma_plot3_red<-
plot(impact3d1_cons_trauma, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Consultas por Causa: Traumatismos")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

cons_resp_plot3_red<-
plot(impact3d_cons_resp, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Consultas por Causa: Sistema Respiratorio")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_trauma_plot3_red<-
plot(impact3d_ratio, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Hospitalizaciones por Causa: Traumatismos, cada 1.000 consultas")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

rate_resp_plot3_red<-
plot(impact3d_ratio_resp, "cumulative")$data %>%
  dplyr::left_join(data15a64_rn[,c("date","rn")], by=c("time"="rn")) %>%
      ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = as.numeric(as.Date("2019-10-18")), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = angle_x, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Hospitalizaciones por Causa: Sistema Respiratorio, cada 1.000 consultas")+
  scale_x_date(date_breaks = "2 weeks",date_labels="%b %d",limits=c(as.Date("2019-08-05"),as.Date("2019-12-23")))

Sys.setlocale(category = "LC_ALL", locale = "spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
plot3b_red<-
plot_grid(
  cons_trauma_plot3_red+scale_y_continuous(breaks=seq(-4000,1500,1000), limits = c(-4000,1500)),
  cons_resp_plot3_red+scale_y_continuous(breaks=seq(-4000,1500,1000), limits = c(-4000,1500)),
  hosp_trauma_plot3_red+scale_y_continuous(breaks=seq(-200,200,100), limits = c(-200,200)),
  hosp_resp_plot3_red+scale_y_continuous(breaks=seq(-200,200,100), limits = c(-200,200)),
  rate_trauma_plot3_red+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3_red+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25), scale = c(.95, .95, .95,.95,.95,.95))+
 draw_label("Diferencias semanales acumuladas en consultas/hospitalizaciones",x=0.005, y=.5, angle=90, size = 17)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

ggdraw(plot3b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))
Figura 2. Diferencias acumuladas entre resultados predichos y observados, en las 10 semanas anteriores y posteriores a la exposición

Figura 2. Diferencias acumuladas entre resultados predichos y observados, en las 10 semanas anteriores y posteriores a la exposición

  jpeg("_figs/_Fig2_final_ESP.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot3b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
## png 
##   2
cairo_ps("_figs/_Fig2_final_ESP.eps", 
           height = 14, width = 23,
           fallback_resolution = 600,
           family= "Times")
           #horizontal=T,
           #paper= "letter", 
           #pagecentre=F)
      ggdraw(plot3b_red)+
         theme(plot.background = element_rect(fill=NA, color = NA))
dev.off()
## png 
##   2
 ggsave("_figs/_FigS2_final_ESP.ps", plot =
          ggdraw(plot3b_red)+
             theme(plot.background = element_rect(fill=NA, color = NA)),
      path = NULL,
       scale = 1, width = 23, height = 14, units = "in",
       device=cairo_ps, 
       fallback_resolution = 600,
       limitsize = F)
vector_bsts_models<-c("impact3d1_cons_trauma", "impact3d_cons_resp","impact3d_hosp_trauma", "impact3d1_hosp_resp", "impact3d_ratio", "impact3d_ratio_resp")
names_outcomes<-c("Consultas por Causa: Traumatismos","Consultas por Causa: Sistema Respiratorio","Hospitalizaciones por Causa: Traumatismos","Hospitalizaciones por Causa: Sistema Respiratorio","Hospitalizaciones por Causa: Traumatismos, cada 1.000 consultas","Hospitalizaciones por Causa: Sistema Respiratorio, cada 1.000 consultas")
df_results_bsts_esp<-data.frame()
for (i in 1:6) {
  x<-vector_bsts_models[i]
dt_bsts<-cbind(
  outcome=names_outcomes[i],
  AE=sprintf("%4.1f",round(get(x)$summary$AbsEffect[1],2)),
  IC95_AE=paste0(sprintf("%4.1f",round(get(x)$summary$AbsEffect.lower[1],2)),", ",sprintf("%4.1f",round(get(x)$summary$AbsEffect.upper[1],2))),
  p=sprintf("%5.3f",round(get(x)$summary$p[1],5)),
  RE=sprintf("%4.1f",round(get(x)$summary$RelEffect[1]*100,2)),
  IC95_RE=paste0(sprintf("%4.1f",round(get(x)$summary$RelEffect.lower[1]*100,2)),", ",sprintf("%4.1f",round(get(x)$summary$RelEffect.upper[1]*100,2)))
  )
  df_results_bsts_esp<-rbind.data.frame(df_results_bsts_esp,dt_bsts)
}
df_results_bsts_esp[which(df_results_bsts_esp$p=="0.000"),"p"]<-"<0.001"

df_results_bsts_esp %>% 
  dplyr::select(-p) %>% 
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 2. Impacto estimado del estallido social de Octubre de 2019 en consultas y hospitalizaciones por causas Traumatismos y Sistema Respiratorio en servicio de Urgencias"),
               col.names = c("","Impacto Absoluto<sup>a</sup>","Intervalo de Credibilidad (95%)", "Impacto Relativo(%)","Intervalo de Credibilidad (95%)"),
               align =c('l',rep('c', 101)),
               escape = F) %>%
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size= 11) %>% 
  kableExtra::add_footnote(c("Estimado como el promedio de la diferencia entre las consultas u hospitalizaciones esperadas y las observadas en el tiempo de exposición"),
                             notation = "alphabet")%>%
  kableExtra::scroll_box(width = "100%", height = "375px")